#include <iostream>
#include <vector>
#include <functional>
#include <cmath>
#include <blitz/array.h>
#include <clay/linear_algebra.h>
#include <system/task_timer.h>
#include <future>

using namespace std;

static constexpr int N = 400;

template <class expr_T>
void evaluate( expr_T&& expr )
{
	for( int i = 0; i < N*N*N; i++)
	{
		//if( i < 100 ) cout << *expr << " "; else
			double x = *expr;
		expr.advance();
		
	}
	cout << endl;
	/*
	cout << *expr << endl;
	expr.advance();
	cout << *expr << endl;
	expr.advance();
	cout << *expr << endl;
	*/
}

template <class T>
T parallel_sum(T* data, int size)
{
  T sum = 0;
  if ( size < 100000 )
    for ( int i = 0; i < size; ++i )
      sum += data[i];
  else {
    auto handle = std::async( parallel_sum<T>, data+size/2, size-size/2);
    sum += parallel_sum(data, size/2);
    sum += handle.get();
  }
  return sum;
}

template<class T, class V> 
struct Accum  
{	// simple accumulator function object
	T* b;
	T* e;
	V val;
	Accum(T* bb, T* ee, const V& v) : b{bb}, e{ee}, val{v} {}
	V operator() () { return std::accumulate(b,e,val); }
};

template <class T>
T comp( vector<T>& v )
{
	if (v.size()<10000) return std::accumulate(v.begin(),v.end(),0.0);

	auto f0 = std::async( Accum<T, T>{&v[0],&v[v.size()/4],0.0} );
	auto f1 = std::async( Accum<T, T>{&v[v.size()/4],&v[v.size()/2],0.0} );
	auto f2 = std::async( Accum<T, T>{&v[v.size()/2],&v[v.size()*3/4],0.0} );
	auto f3 = std::async( Accum<T, T>{&v[v.size()*3/4],&v[v.size()],0.0} );

	return f0.get()+f1.get()+f2.get()+f3.get();
}
	
int main(int argc, char** argv)
{
	vector<double> v = {1.0, 3.4, 5.6};
	array<complex<double>, 1> arr{3};
	arr = { {1,2}, {3,4}, {5,6} };
	array<double, 1> arr2{3};
	arr2 = { 7, 8, 9 };
	complex<double> con_v = {1.0, 2.0};
	
	//expression<decltype(arr2 + arr2)> expr_arr(arr2 + arr2);
	//expression_func<decltype(expr_arr), decltype(func)> expr_func(expr_arr, func);
	//expression<decltype(expr_func)> expr_end(expr_func);
	
	//auto sin_expr = sin(arr2 + arr2);
	//cout << con_v + sin(con_v+con_v) + sin(2.0) << endl;
	//evaluate( arr + sin(arr+arr) + sin(2.0));
	
	//evaluate( arr[0] / arr2 + 4.3 / arr2 * arr[1] - 0.3*arr2);

//	array<tiny_vector<double,3>, 1> vecarr{3};
//	vecarr = { {0,0,0}, {1,1,1}, {2,2,2} };
	auto func = []( double a, double b, double c ) { return a* b + c; };
	//auto expr = create_expression_func( func, arr2, 2.0, 1.0 );
	//evaluate( expr );
	auto axpy = create_expr_func_adapter(func);
	//evaluate( func_adp_multi(arr2, 2.0, 1.0) );
	//evaluate( spx::sin2(arr2) );
	//evaluate( spx::pow( spx::pow( arr2, 2.0 ), 2.0) );
	//evaluate( pow( arr2, 2.0 ) * pow( arr2, 2.0 ) );
	//auto func_adp = declare_expr_func_adapter(func);
//	evaluate( func_adp(vecarr) );
	//evaluate( spx::polar( arr2*2.0, arr2 ) );
		
	array<double, 3> arr3{300};
	arr3 = 2.0;
	tiny_vector<int, 3> a = {1,2,3};
	auto b = a + 4;
	cout << b << endl;
	auto task = [&] {
		//arr3 += spx::pow(arr3, arr3) + 2.0; //spx::pow(arr3, 2.0); 
		//evaluate( spx::pow(arr3, 3.0) * 2.0 );
		//evaluate( spx::log(arr3) + 2.0 );
		//evaluate( axpy(2.0, arr3, 1.0) );
		//evaluate( arr3 < 10.0 );
		//arr3 += arr3 < 10.0;
		//evaluate( -arr3 + arr3 + 4.0 + spx::sin(arr3) );
		//arr3 = -arr3 + arr3 + 4.0 + spx::pow(arr3, 2.1);
		array<double, 3> arr4( -arr3 + arr3 + 4.0 + spx::sin(arr3) );
	};
	task_timer<> timer;
	timer( task );
	cout << timer.get_total_time() << endl;
	cout << arr3(90,90,90) << endl;
	
	BZ_USING_NAMESPACE(blitz)
	timer.set_zero();
	Array<double, 3> Arr3(300, 300, 300);
	Arr3 = 2.0;
	//Array<double, 3> Arr4 = Arr3 + 2.0;
	auto task2 = [&] {
		//Arr3 += pow(Arr3, Arr3) + 2.0; //pow(Arr3, 2.0); 
		//evaluate( pow(Arr3, 3.0) + 2.0 );
		//evaluate( log(Arr3) + 2.0 );
		//evaluate( 2.0 * Arr3 + 1.0 );
		//evaluate( Arr3 < 10.0 );
		//Arr3 += Arr3 < 10.0;
		//evaluate( -Arr3 );
		//evaluate( Arr4 + 1.0 );
		//evaluate( -Arr3 + Arr3 + 4.0 + sin(Arr3) );
		//Arr3 = -Arr3 + Arr3 + 4.0 + pow(Arr3, 2.1);
		Array<double, 3> Arr4( -Arr3 + Arr3 + 4.0 + sin(Arr3) );
	};
	timer( task2 );
	cout << timer.get_total_time() << endl;
	cout << Arr3(90,90,90) << endl;
	
	/*
	int size = 100000000;
	vector<double> a(size);
	double sum = 0.0;
	for( auto& v : a ) v = 1.0;
	task_timer<> timer2, timer4;
	auto task3 = [&sum, &a, size] () { 
		for(int i =0; i< 50;i++)
			//sum = parallel_sum( &a[0], size ); 
			sum = comp( a );
	};
	timer2( task3 );
	cout << sum << endl << timer2.get_total_time() << endl;
	
	auto task4 = [&sum, &a] () {
		
		#pragma omp parallel for
		for(int i =0; i< 50;i++)
		{
			double sum_ = 0.0;
			for( int n = 0; n < a.size(); n++ ) sum_ += a[n]; 
			sum = sum_;
		}
	};
	timer4( task4 );
	cout << sum << endl << timer4.get_total_time() << endl;
	*/
}
